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ABSTRACT 


The  method  of  high-resolution  spectral  analysis  has  been 
applied  to  one-  and  two-dimensional  wavenumber  spectra.  In 
this  method  a  prediction  operator  is  used  to  extend  the  length 
of  the  spatial  correlation  functions  in  order  to  obtain  higher 
resolution  of  the  resulting  spectra  in  wavenumber  space.  It 
is  most  useful  in  the  case  of  small  arrays,  where  the  spatial 
correlation  functions  are  too  short  to  give  adequate  resolu¬ 
tion  for  the  ordinary  frequency-wavenumber  spectra.  A  theo¬ 
retical  description  of  the  method,  its  procedures,  and  its 
relation  to  ordinary  wavenumber  spectra  is  given.  In  partic¬ 
ular,  an  analytical  expression  is  derived  for  the  resolution 
and  the  response  of  the  process  in  the  case  of  an  input  plane 
wave,  where  the  high-resolution  method  is  averaged  over  all 
reference  sensors. 

Disadvantages  of  the  high-resolution  method  are  discussed, 
in  particular  an  unrecoverable  distortion  of  the  true  amplitude 
spectrum.  Finally,  examples  are  given  for  data  recorded  at  the 
Wichita  Mountains  Seismic  Observatory  surface  array  and  the 
Apache,  Oklahoma,  vertical  array. 
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INTRODUCTION 


Burg  (1967)  has  shown  that  high-resolution  power  spectra 
can  be  computed  from  relatively  short  correlation  functions  by 
first  least-squares  predicting  larger  correlation  lags  from 
those  already  known.  The  resulting  high-resolution  power  spec¬ 
trum  is  inversely  proportional  to  the  power  spectrum  of  the  cor¬ 
responding  prediction  error  operator.  However,  with  the  advent 
of  the  Cooley-Tukey  fast  Fourier  transform  algorithm  and  related 
methods  for  estimating  power  spectra,  the  high-resolution  tech¬ 
nique  for  ordinary  power  spectra  is  no  longer  economical.  It 
is  now  feasible  to  compute  power  spectra  as  if  they  were  Fourier 
transforms  of  the  100%- lag  correlation  functions;  various  kinds 
of  smoothing  can  then  be  applied  to  reduce  the  spectra  to  the 
desired  length  and  stability. 

We  feel  that  the  most  promising  application  of  the  high- 
resolution  procedure  lies  in  estimating  frequency-wavenumber 
spectra.  Here  the  length  of  the  spatial  correlation  functions, 
which  ultimately  determines  the  resolution,  is  limited  by  the 
size  of  the  array.  The  high  resolution  method  is  therefore 
the  appropriate  way  to  obtain  frequency-wavenumber  spectra 
which  are  sharp  in  wavenumber  space  from  data  recorded  over 
small  arrays. 

Ordinary  frequency-wavenumber  spectra  (Spieker  et  al. , 

1961)  are  computed  directly  from  the  spectral  matrix  elements 
by  the  expression: 

N  N  r  "1 

P(w,k)  =  l  l  S..(w)  exp  -ik.(x.  -  x . ) 

i=l  j=l  1J  L  JJ 

In  the  case  of  two-dimensional  arrays  (for  example,  LASA)  P  is 
a  function  of  three  variables:  w,  kx,  and  k^ .  The  standard  way 
of  displaying  this  function  is  to  produce  a  succession  of  con¬ 
tour  plots,  each  at  a  particular  frequency,  with  P  contoured 

against  k  and  k  .  On  the  other  hand,  for  one-dimensional 
~  x  ~y 


l 


arrays  (for  example,  arrays  of  instruments  in  deep  wells)  only 
the  vertical  wavsnumber  is  involved  and  P  can  be  contoured 
against  u  and  k2  in  one  plot.  Loci  of  constant  velocity 

co  k 


are  circles  about  the  k*,  ky  origin  as  shown  in  Figure  la  or 
lines  of  constant  slope  as  shown  in  Figure  lb. 

A  measure  of  the  resolution  or  capability  of  the  array  is 
its  response  to  an  infinite-velocity  plane  wave.  If  the  instru¬ 
ment  frequency  response  is  neglected,  the  spectral  matrix  of 
the  wave  is  constant.  Therefore  this  function,  usually  called 
the  array  response  function  (Spieker  et  al. ,  1961), i»  given  by: 

N  N 


R(*>  ~  l  l  exp  f-ik*(x.  -  x.)l 

i=l  j=l  L  ~  ~i  J 


Any  single  f inite-velocity  plane  wave  of  the  form 

exp  [  i  ( kQ  •  x  -  cot )  3 

will  produce  a  frequency-wavenumber  spectrum  consisting  of  the 
array  response  function  centered  on  the  value  k  modulated  in 
frequency  by  the  corresponding  frequency  power  spectrum. 

The  high-resolution  technique  consists  of  using  the  spec¬ 
tral  matrix  of  the  data  to  design  a  multichannel  prediction 
error  operator.  High-resolution  frequency-wavenumber  spectra 
are  then  computed  from  the  spectral  matrix  of  this  multichan¬ 
nel  filter  (Haney,  1967).  The  appropriate  matrix  equation  is: 

tS(co)  +  cl  ]  £m(w)  =  gm 

Here  S<w>  is  the  spectral  matrix  of  the  data,  I  is  an  N-by-N 
identity  matrix,  £m(w)  is  the  filter  vector,  and  gm  is  the 
nghthand  side  vector  consisting  of  N-l  zeros  and~a  one  in 
the  m-  row.  This  equation  defines  the  filter  which  is  opti¬ 
mum,  in  the  least-squares  sense,  for  "whitening"  the  m^ 
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channel.  Another  way  of  interpreting  equation  (5)  is  to  rec¬ 
ognize  it  as  a  multichannel  Wiener- Hopf  equation  with  S(w)  the 
noise  spectral  matrix,  cl  the  signal  model,  and  g  proportional 
to  the  cross  spectra  between  the  desired  output,  white  noise,  and 
the  input.  The  constant  c  is  a  convenient  means  of  adjusting  the 
signal-to-noise  ratio  of  the  input. 

The  spectral  matrix  of  the  filter  is  then 

Hm(w)  =  E  [fm(w)fm(w)+]  (6) 

where  the  dagger  indicates  Hermitian  conjugation  and  the  brackets 
specify  that  expected  values  are  to  be  taken.  High-resolution 
frequency-wavenumber  spectra  of  the  original  data  channels  are 
inversely  proportional  to  the  ordinary  frequency-wavenumber 
spectra  of  this  spectral  matrix: 

P^(w,k)  =1^  H^.(w)exp  (-i‘-j )  }  (7) 

In  practice,  we  have  found  that  high-resolution  frequency- 
wavenumber  spectra  computed  this  way  (from  single-output  whiten¬ 
ing"  filters)  sometimes  suffer  an  angular  distortion  (Lintz,  1958). 
In  other  words,  the  resolution  in  wavenumber  space  is  a  function 

of  k-k  .  Furthermore,  it  is  undesirable  to  have  the  technique 
~  ~o 

depend  on  the  choice  of  the  reference  data  channel.  This  diffi¬ 
culty  can  be  circumvented  by  using  the  average  spectral  matrix 
of  the  filters,  H,  computed  by  averaging  over  all  reference  datS 
channels.  This  procedure  is  also  a  convenient  means  for  taking 
the  expected  value. 

Thus : 


H(w) 


j  Hm(w)  =  l  fm(u>)fm(u>)  + 
m=l  m=l 


(8) 


Since  each  filter  vector  fm(w)  is  merely  the  m —  column  of  the 
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inverse  matrix: 

[S(w)  +  cl]  ( 

the  required  averaged  spectral  matrix  HCw)  is  easily  shown  to  be 


(10) 

This  modification  of  the  original  procedure  is  used  for  the 
examples  given  below.  It  should  be  noted  that  computing  high- 
resolution  frequency-wavenumber  spectra  this  way  is  only  negli¬ 
gibly  more  expensive  than  computing  ordinary  frequency-wavenumber 
spectra.  Both  techniques  require  computing  the  spectral  matrix. 

The  only  other  computations  required  by  the  high-resolution  tech¬ 
nique  are  the  addition  of  some  white  noise,  one  Hermitian  matrix 
inversion,  and  one  Hermitian  matrix  multiplication. 

The  actual  improvement  in  resolution  gained  by  using  the 
averaged  high-resolution  technique  can  be  demonstrated  by  applying 

it  to  a  plane  wave.  Neglecting  the  frequency  part,  the  spectral 
matrix  of  this  wave  is: 


sij(u)  *  exP  [i5o*<*i-*j]| 

which  can  be  factored  into  the  outer  product  of  two  vectors : 

S  (u> )  =  q  q  + 

where  the  vector  q  is  given  by: 


|9_|i  =  exP  C ikQ •  xi ) 
The  matrix  to  be  inverted  is: 


(  q  q+  +  cl) 


which  is  equal  to: 


— —  (  T  H 

c  V1  c  +  N 


Squaring  this  matrix  gives,  to  terms  of  order  c3, 


I  -  2S  x  NS 

C+N  7^7 


(id 


(12) 


(13) 


(14) 


(15) 


*  E's  -  -S-)+ 


0(c)  (16) 
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Substituting  equation  (16)  into  equation  (7)  gives  the  result: 


P*  (k)  - - 

N  -  P(k) 


(17) 


2 

Now  when  POO  reaches  a  maximum  value  of  N  ,  P'(k)  has  the  value 
N. 

However  when  POO  has  fallen  off  to  N2/a,  P'OO  is: 


P'  (k) 


(18) 


which,  by  choosing  c  to  be  much  less  than  N; 

c  <<  N  (19) 

makes  P5(k)  very  small  indeed.  The  resolution  of  this  technique 
therefore  is  dependent  upon  both  N,  the  number  of  data  channels, 
and  the  constant  c. 
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PROCEDURES 


All  of  the  spectral  matrix  elements  computed  for  the  examples 
given  below  were  done  so  by  the  Cooley-Tukey  method  (McCowan, 1967 ) . 
The  data  was  Fourier  transformed  togetl  er  with  a  sequence  of  zeros 
equal  to  its  entire  length  and  then  cross-multiplied  to  produce 
the  raw  auto-  and  cross-spectra.  These  raw  spectra  were  then  sub¬ 
jected  to  a  number  of  smoothing  operations.  Each  operation  con¬ 
sisted  of  convolving  the  spectra  with  the  three-point  operator 
(1/4,  1/2,  1/4)  and  taking  every  other  point  of  the  result.  The 
effects  of  this  window  are  described  elsewhere  (McCowan,  1967). 
Typically  the  original  data  channels  were  of  the  order  of  4096 
digital  points  for  noise  samples  and  256  digital  points  for  sig¬ 
nals,  both  sampled  at  twenty  points  per  second.  This  gave  4097- 
point  and  257-point  raw  spectra,  respectively,  which  were  then 
smoothed  and  reduced  to  65  and  33  actual  frequency  values  respect¬ 
ively.  Thus  the  number  of  degrees  of  freedom  included  in  each 
spectral  estimate  was  64  in  the  case  of  noise  samples  and  8  in  the 
case  of  signals. 

Both  to  alleviate  the  problem  caused  by  instrument  gain  ano¬ 
malies  (Baldwin,  1964)  and  the  question  of  physical  units  in 
choosing  the  constant  c  in  equation  (5),  we  have  used  normalized 


spectral  matrices  (coherence  matrices)  throughout.  Our  spectral 
matrix  elements  are  thus: 

E{F . (w)F. (w)*} 

S..U)  =  - - - — 2 - 5 -  (20) 

3  E{  |Fi(w)  |  }  E{  |  Fj  (w)  p} 


Otherwise,  choosing  the  constant  c,  and  hence  determining  the 
resolution  of  the  technique,  would  depend  on  the  actual  magnitude 
of  the  auto- spectra . 

In  the  examples  below  we  have  chosen  frequency  in  cps  instead 
of  angular  frequency  as  an  independent  variable.  Consequently, 
we  have  re-defined  the  wavenumber  k  as: 

|k|  =  -j —  (21) 
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so  that  equation  (2)  now  reads: 


This  was  done  purely  as  a  matter  of  convenience  in  reading  the 
countour  plots. 

Finally,  the  contour  plots  displayed  below  were  computed  on 
a  grid  of  21- by- 21  k  and  k  values  for  the  surface  array  and 
21-by-33  k  and  frequency  values  for  the  vertical  array.  Both 
were  then  linear-interpolated  in  each  direction,  to  a  lOO-by-60 
grid  to  be  printed  by  the  computer  printer.  These  pages  were  then 
traced  to  produce  the  figures.  Any  details  smaller  than  1<  20  of 
the  whole  figure  are  therefore  strongly  colored  by  the  linear 
interpolation. 


EXAMPLES 


Two  events  (P  wave  signals  summarized  in  Table  1)  and  samples 
of  noise  immediately  preceding  them  were  selected  as  examples  of 
high-resolution  frequency-wavenumber  data  processing.  The  first 
event  and  its  corresponding  noise  sample  were  used  to  illustrate 
the  application  of  the  technique  to  a  two-dimensional  surface 
array,  in  this  case  the  WMSO  array,  while  the  remaining  data  were 
used  to  illustrate  results  from  a  one-dimensional  vertical  array, 
the  APOK  deep  well.  Both  arrays  are  descibed  by  the  coordinates 
of  their  elements,  which  are  listed  in  Table  2,  and  their  array 
response  functions,  which  are  shown  as  contour-plots  in  Figures  3 
and  8  respectively.  A  map  of  the  WMSO  array  is  included  as 
Figure  2. 


Ordinary  frequency-wavenumber  spectra  of  both  pairs  of  noise 
samples  and  signals  are  shown  in  Figures  4,  6,  9a,  and  10a.  Figures 
4  and  9a  are  the  noise  samples  and  Figures  6  and  10a  are  the  sig¬ 
nals.  The  WMSO  noise  seems  to  be  coming  from  the  northeast  with  a 
velocity  ranging  between  3  and  6  km/sec.  This  characteristic  behav¬ 
ior  has  been  noted  by  other  investigators  (Rizvi  et  al.,  1967). 

On  the  other  hand,  the  WMSO  signal  appears  to  be  coming  from  the 
souch  or  southeast  with  a  velocity  somewhere  between  16  and  20 
km/sec.  The  APOK  noise  sample  displays  a  tilted  pattern  also 
characteristic  of  the  site.  One  explanation  for  this  phenomenon 
is  suggested  by  the  geological  structure  in  the  vicinty  of  the 
well  (Geotechnical  Corp.,  1964),  which  consists  of  Paleozoic  in- 


trusives  and  sediments  dipping  30-40°  towards  the  northeast.  The 
APOK  signal  is  very  strange  indeed,  composed  of  a  wave  travelling 
down  the  well  with  an  apparent  vertical  velocity  of  approximately 
3.5  km/sec  and  a  horizontally  travelling  wave  with  infinite  ver¬ 
tical  velocity.  One  possible  explanation  of  this  result  is  again 
suggested  by  the  dip  of  the  structure  (Sax,  1967).  The  incoming  P 
wave  might  be  refracted  by  the  tilted  interface  so  that  the  wave 
front  would  be  a  leaking  suface  wave  when  it  reached  the  well 


High-resolution  frequency-wavenumber  spectra  are  shown  for 
all  these  data  samples  in  Figures  5,  7,  9b,  and  10b.  These 
were  computed  with  a  signal-to-noise  ratio  of  which  seemed 
to  be  a  good  balance  between  high  resolution  and  excessive  dis¬ 
tortion.  For  the  WMSO  examples,  the  technique  decreased  the 
width  of  the  -3db  contour  by  more  than  a  factor  of  two  for  both 
the  cases  of  noise  and  signal.  The  results  are  more  difficult 
to  assess  for  the  vertical  array  data  because  of  the  whitening 
in  frequency  done  by  the  filters.  However,  for  the  noise  sample 
at  0.25  cps,  a  decrease  by  a  factor  of  two  in  the  width  of  the 
-3  db  contour  can  again  be  seen. 

A  marked  disadvantage  and  a  difficulty  still  to  be  solved 
in  the  application  of  this  technique  is  illustrated  in  Figure  10b. 
The  high-resolution  technique,  when  applied  to  short  samples  of 
data,  can  produce  spurious  peaks  in  the  frequency-wavenumber 
spectra.  These  are  due  to  the  limited  smoothing  done  when  com¬ 
puting  spectra  from  short  samples  of  data.  Figure  10b  shows  sev¬ 
eral  of  these  peaks.  This  is  regarded  as  an  unrecoverable  dis¬ 
tortion  and  in  a  certain  sense,  an  inevitable  result  of  using  a 
high-gain  procedure. 
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Figure  8.  APOK  array  response 
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Table  1 


EVENT  DATA 

WMSO  Event 

APOK  Event 

Latitude 

41.4°  S. 

52.0°  N. 

Longitude 

88.6°  W. 

175.0°  W. 

Origin  Time 

14: 51:55.0 

19:05:04.1 

Depth  (km) 

33 

67 

Magnitude 

4.9 

5.2 

Region 

W.  Chile  Rise 

Aleutians 

Data  taken  from  USCGS  Preliminary  Determination  of  Epicenters 


Table  2 


SENSOR  COORDINATES 

WMSO  SURFACE  ARRAY 


Sensor 

X(km. ) 

Y(kin.  ) 

V1 

0.00 

0.00 

V  2 

0.67 

0.74 

V  3 

0.94 

-0.17 

V4 

0.27 

-0.70 

V5 

-0.99 

0.12 

V6 

-0.34 

0.71 

V7 

0.08 

1.41 

V8 

1.35 

1.33 

V9 

1.68 

0.34 

< 

H 

O 

1.19 

-1.02 

V11 

-0.54 

00 

r- 

o 

V12 

-1.  51 

-0.59 

V13 

-1.29 

0.95 

APOK  VERTICAL  ARRAY 

Sensor _  Depth  (km.) 


DW5 

0.02 

1.68 

dw3 

1.99 

dw2 

2.29 

DU1 

2.91 

Security  ClMEification 


DOCUMENT  CONTROL  DATA  •  RAD 

(Saeurtly  elMilHcalltn  of  Hilt,  body  of  afeofroef  and  Indaalng  annotation  mual  bo  antarad  mfion  lha  aaarall  tapoit  it  a  lata  Iliad) 


I.  ORIGIN  ATIN  0  ACTIVITY  (Cotparala  aulhar)  la.  RIRORT  (CCURITV  C  LAtdRIC A f  |4U 

TELEDYNE,  INC.  Unclassified 

ALEXANDRIA,  VIRGINIA  22314  77 onou'p - 


t  REPORT  TITLE 


HIGH-RESOLUTION  FREQUENCY-WAVENUMBER  SPECTRA 


J  i  MJLiiri  f'Mi  iWki 


Scientific 


|.  AUTHOHW  (L»ml  namt,  tint  Italia,  flnfiimtj 


Lintz,  P.R.,  McCowan,  D.W. 


I  REPO  NT  OATE 

16  February  1968 


•  a.  CONTRACT  OR  «RANT  NO. 

F  33657-67-C-1313 

k  RROJCCT  NO. 

VELA  T/6702 


ARPA  Order  624 


la.  TOTAL  NO.  or  PAIN  76.  NO.  OR  RKRI 

26  9 


to.  ORIOINATOR't  RIRORT  NUMWBRfS) 

SDL  Report.  No.  206 


*6  W%.nr,r0nr  MOW  <A  "r  OU,m'  niJ<nk<r*  ba  aaal0tad 


11  SUPPLEMENTARY  NOTES 


10  AVAIL  ABILITY/LIMITATION  NOTICES 

This  document  is  subject  to  special  export  controls  and  each  trans 
mittal  to  foreign  governments  or  foreign  national  may  be  made  only 
with  prior  approval  of  Chief,  AFTAC. _ 


1*.  SPONSORINO  MILITARY  ACTIVITY 

ADVANCED  RESEARCH  PROJECTS  AGENCY 
NUCLEAR  TEST  DETECTION  OFFICE 
WASHINGTON ,  D.C. 


is  abstract  The  method  of  high-resolution  spectral  analysis  has  been  ap¬ 
plied  to  one-  and  two-dimensional  wavenumber  spectra.  In  this  meth¬ 
od  a  prediction  operator  is  used  to  extend  the  length  of  the  spatial 
correlation  functions  in  order  to  obtain  higher  resolution  of  the 
resulting  spectra  in  wavenumber  space.  It  is  most  useful  in  the 
case  of  small  arrays,  where  the  spatial  correlation  functions  are 
too  short  to  give  adequate  resolution  for  the  ordinary  frequency- 
wavenumber  spectra.  A  theoretical  description  of  the  method,  its 
procedures,  and  its  relation  to  ordinary  wavenumber  spectra  is  given. 
In  particular,  an  analytical  expression  is  derived  for  the  resolu¬ 
tion  and  the  response  of  the  process  in  the  case  of  an  input  plane 
wave,  where  the  high-resolution  method  is  averaged  over  all  ref¬ 
erence  sensors . 

Disadvantages  of  the  high-resolution  method  are  discussed,  in 
particular  an  unrecoverable  distortion  of  the  true  amplitude  spec¬ 
trum.  Finally,  examples  are  given  for  data  recorded  at  the  Wichita 
Mountains  Seismic  Observatory  surface  array  and  the  Apache,  Oklahoma, 
vertical  array.  ' 
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